/****************************************************************
*****************************************************************/

*by Xiaodong Fan, fanxiaodong@gmail.com

#delimit ;

cap log close;
clear all;
drop _all;
set more 1;
pause on;

log using log_MSM_elas.log, replace;

global gpicsdir "pics_elas";
global gdtadir ".";

include do_globals.do;

if ( (strpos("`c(pwd)'", "8") > 0) | (strpos("`c(pwd)'", "9") > 0) ) {;
   global gvpos = 11;
};
else {;
   global gvpos = 11;
};

***********************************************************;
* first of all, parameter estimates and standard errors;
***********************************************************;
*import delimited pars_mapping.csv, clear;
insheet using pars_mapping.csv, comma names clear;
keep if [_n]<=${gvparnum};
sort incode;
save dta_pars_mapping.dta, replace;

infile par using "txt_parin9.txt", clear;
keep if [_n]<=${gvparnum};
destring, replace;
gen incode = [_n];

sort incode;
merge 1:1 incode using dta_pars_mapping.dta, nogenerate;
sort incode;
save dta_10_pars.dta, replace;

* standard errors;
local lvparse "N2";
infile incode var se using "txt_par_se_`lvparse'.txt", clear;
keep if [_n]<=${gvparnum};
destring, replace;
compress;
sort incode;
merge 1:1 incode using dta_10_pars.dta, nogenerate;
sort incode;

sort incode;
save dta_10_pars.dta, replace;

** the consumption shifter parameters have already been scaled by 10^i in the Fortran program;
if (1<0) {;
    prog_parSSAx1000 par se;
    sort inpaper;
    export excel incode inpaper par se using "table4paper${alt_name}", sheet("pars_se${alt_name}") sheetreplace firstrow(variables);
};

/***********************************************
  Moments generated in Fortran program
************************************************/
local lvmoments lfpr lnw sdlnw lnw_fd alfpr1to0 alfpr0to1 C ssa lfpr_diff_E2G lfpr_diff_G2B lfpr_diff_B2D lfpr_pt plnw41;
local lvAH0 A H C I labor inc V0 V1 probV0;
local lvAH A H C_profile I labor inc V0 V1 probV0 Vp probVp;

* baseline model;
prog_MSM_moments_infile "txt_MSM_Moments.txt" "";
ren *_f *;

sort t;
replace lnw_fd = . if t<${gvdata0} | t>${gvdata10};
gen byte elast = 0;
sort elast t;
save ${gdtadir}/dta_MSM_elas.dta, replace;

* merge to get other profiles in the baseline model;
use ${gdtadir}/dta_MSM_Profiles.dta, clear;
keep if ihet == 0;
keep t `lvAH0';
ren C C_profile;
gen byte elast = 0;
sort elast t;
merge 1:1 elast t using ${gdtadir}/dta_MSM_elas.dta, update replace;
tab _merge;
drop _merge;
sort elast t;
save ${gdtadir}/dta_MSM_elas.dta, replace;

* merge with the elas file;
prog_MSM_moments_infile "txt_MSM_Moments_elas.txt" "elast" "`lvAH'";
ren *_f *;

sort t;
replace lnw_fd = . if t<${gvdata0} | t>${gvdata10};
append using  ${gdtadir}/dta_MSM_elas.dta;
compress;

sort t elast;

* plot graph;
preserve;
   local lvx1 lfpr lnw lnw_fd `lvAH';
   keep t elast `lvx1';
   reshape wide `lvx1', i(t) j(elast);
   sort t;
   foreach iv in lfpr H I lnw {;
      local lvyscale "";
      local lvylabel "";
      local lvpos = 12;
      if ("`iv'"=="lfpr") {;
         local lvytitle "Labor Force Participation Rate";
      };
      else if ("`iv'"=="lnw") {;
         local lvytitle "Log Wages";
      };
      else if ("`iv'"=="lnw_fd") {;
         local lvytitle "Log Wages (FD)";
      };
      else if ("`iv'"=="C") {;
         local lvytitle "Consumption";
      };
      else if ("`iv'"=="I") {;
         local lvytitle "Investment";
         local lvpos = 5;
         local lvyscale " yscale(range(-0.01 0.005))";
         local lvylabel "ylabel(-0.01(0.005)0.005)";
      };
      else if ("`iv'"=="H") {;
         local lvytitle "Human Capital";
      };
      else {;
         local lvytitle "`iv'";
      };

      foreach iy in 25 35 40 60 65 {;
         gen d_`iv'`iy' = `iv'`iy' - `iv'0;
      };

      line d_`iv'25 d_`iv'40 d_`iv'60 t if t<=${gvt_elas}, lpattern("l" "_" "-" "-.") `lvyscale' `lvylabel'
          xtitle("Age") ytitle("`lvytitle'") graphregion(color(white)) xlabel(20(10)${gvt_elas} ${gvt_elas},grid) 
          legend(pos(`lvpos') ring(0) col(1) symxsize(*1) order(1 "Shock at 25" 2 "Shock at 40" 3 "Shock at 60"));
      graph export ${gpicsdir}/eps_elas_shock_`iv'.eps, replace;
 
      line d_`iv'25 d_`iv'35 d_`iv'60 d_`iv'65 t if t<=${gvt_elas}, lpattern("l" "_" "-" "-.") `lvyscale' `lvylabel'
          xtitle("Age") ytitle("`lvytitle'") graphregion(color(white)) xlabel(20(10)${gvt_elas} ${gvt_elas},grid) 
          legend(pos(`lvpos') ring(0) col(1) symxsize(*1) order(1 "Shock at 25" 2 "Shock at 35" 3 "Shock at 60" 4 "Shock at 65"));
      graph export ${gpicsdir}/eps_elas_shock4_`iv'.eps, replace;
 
   };
restore;

foreach iv in `lvmoments' `lvAH' {;
   sort t elast;
   by t: gen `iv'_b = `iv'[1] if [_n]>1;
};
drop if elast == 0;
compress;
save ${gdtadir}/dta_MSM_elas.dta, replace;


/*********************************************************
* calculating different components of elasticity;
*********************************************************/
* get the parameters;
* 1 Leisure: Standard Deviation of Shock;
* 6 Leisure: Mean of Intercept;
preserve;
   use dta_10_pars.dta, clear;
   keep if incode==6 | incode==1;
   keep incode par;
   sum par if incode==6;
   global gva0 = r(mean);
   sum par if incode==1;
   global gvaepsilon = r(mean);
restore;

/*****************/
preserve;
   keep if t==elast;
   gen V01d = V0 - V1;
   gen V01d_b = V0_b - V1_b;
   sum V01d V01d_b;
   *pause;
   gen lnV01d = log(V01d);
   gen lnV01d_b = log(V01d_b);
   gen lnphi_lnV01d = log( normal((lnV01d-${gva0})/${gvaepsilon}) );
   gen lnphi_lnV01d_b = log( normal((lnV01d_b-${gva0})/${gvaepsilon}) );

   gen elas1 = (lnphi_lnV01d - lnphi_lnV01d_b) / (lnV01d - lnV01d_b);
   gen elas2 = (lnV01d - lnV01d_b) / log(1.1);
   gen elas = elas1 * elas2;
   keep t elas elas1 elas2 *V01d*;
   sort t;
   compress;
   save ${gdtadir}/dta_MSM_elas_calculated.dta, replace;

   * plot graph;
   sort t;
   line elas elas1 elas2 t if t<=${gvt_elas}, lpattern("l" "_" "-" "-.")  xtitle("Age") ytitle("Elasticities") graphregion(color(white))
          xlabel(20(10)${gvt_elas} ${gvt_elas},grid)
          legend(pos(11) ring(0) col(1) colgap(*1) symxsize(*1) order(1 "Elasticity" 2 "dln(Phi)/dlnv" 3 "dlnv/dlnw"));
   graph export ${gpicsdir}/eps_elas_parts.eps, replace;

   line elas1 elas2 t if t<=${gvt_elas}, lpattern("l" "_" "-" "-.")  xtitle("Age") ytitle("Elasticities") graphregion(color(white))  ||
        line lnV01d t if t<=${gvt_elas}, lpattern("..") yaxis(2) xlabel(20(10)${gvt_elas} ${gvt_elas},grid)
          legend(pos(12) ring(0) col(1) colgap(*1) symxsize(*1) order(1 "dln(Phi)/dlnv" 2 "dlnv/dlnw" 3 "lnV01d"));
   graph export ${gpicsdir}/eps_elas_parts_lnV01d.eps, replace;

   line elas1 elas2 t if t<=${gvt_elas}, lpattern("l" "_" "-" "-.")  xtitle("Age") ytitle("Elasticities") graphregion(color(white))  ||
        line lnphi_lnV01d t if t<=${gvt_elas}, lpattern("..") yaxis(2) xlabel(20(10)${gvt_elas} ${gvt_elas},grid)
          legend(pos(9) ring(0) col(1) colgap(*1) symxsize(*1) order(1 "dln(Phi)/dlnv" 2 "dlnv/dlnw" 3 "lnPhi"));
   graph export ${gpicsdir}/eps_elas_parts_lnPhi.eps, replace;

restore;
/*****************/


/*********************************************************
* calculating elasticities;
*********************************************************/

* calculate Marshallian (uncompensated) elasticity;
preserve;
   keep if t==elast;
   gen me0 = (log(lfpr)-log(lfpr_b))/log(1.1);
   gen me2 = (log(lfpr)-log(lfpr_b))/(lnw-lnw_b);
   *pause;
   keep t me0 me2;
   foreach iv in me0 me2 {;
   *   replace `iv' = . if `iv' > 20 | `iv'<-4;
      replace `iv' = . if `iv'<-4;
   }; 

   sort t;
   merge 1:1 t using ${gdtadir}/dta_MSM_elas_calculated.dta;
   tab _merge;
   drop _merge;
   compress;
   save ${gdtadir}/dta_MSM_elas_calculated.dta, replace;
restore;

* calculate IES;
* ies0: denominator is log(1.1)---changes in the Human capital rental rate;
* ies2: denominator is ((lnw2/lnw1)-(lnw_b2/lnw_b1))---changes in wages;
preserve;
   keep if t==elast | t==(elast-1);
   sort elast t;
   by elast: gen tn = [_n];
   by elast: gen tnt = [_N];
   drop if tnt<=1;
   drop tnt;
   keep elast tn lfpr lfpr_b lnw lnw_b;
   reshape wide lfpr lfpr_b lnw lnw_b, i(elast) j(tn);
   gen ies0 = (log(lfpr2/lfpr1)-log(lfpr_b2/lfpr_b1))/log(1.1);
   gen ies2 = (log(lfpr2/lfpr1)-log(lfpr_b2/lfpr_b1))/((lnw2-lnw1)-(lnw_b2-lnw_b1));
   ren elast t;
   keep t ies0 ies2;

   foreach iv in ies0 ies2 {;
   *   replace `iv' = . if `iv' > 20 | `iv'<-4;
      replace `iv' = . if `iv'<-4;
   }; 

   sort t;
   merge 1:1 t using ${gdtadir}/dta_MSM_elas_calculated.dta;
   tab _merge;
   drop _merge;
   sort t;
   compress;
   save ${gdtadir}/dta_MSM_elas_calculated.dta, replace;

   * plot graph of me and ies;
   local lvlastT = 70;
   sum me0 if t<=`lvlastT';
   local lvx0_max = ceil(r(max)*2) * 0.5;
   local lvx0_min = floor(r(min)*10) * 0.1;
   local lvx0_step = `lvx0_max'/5;
   sum me2 if t<=`lvlastT';
   local lvx2_max = ceil(r(max)*2) * 0.5;
   local lvx2_min = floor(r(min)*10) * 0.1;
   local lvx2_step = `lvx2_max'/5;
   local lvelasymax = max(`lvx0_max',`lvx2_max');
   if (${gvpos}==11) {;
      local lvelasymin = 0.2;
   };
   else {;
      local lvelasymin = min(`lvx0_min',`lvx2_min');
   };
   local lvelasstep = `lvelasymax'/5;
   *local lvelasymax = 1.8;

   * ylabel(0.2 0.5(0.5)1.5 1.8);

   sort t;


   line me0 ies0 me2 ies2 t if t<=`lvlastT', lpattern("l" "_" "-" "_.")  xtitle("Age") ytitle("Elasticity") graphregion(color(white))
          xscale(range(18 `lvlastT')) xlabel(20(10)60 `lvlastT',grid) 
          yscale(range(`lvelasymin' `lvelasymax')) ylabel(0(`lvelasstep')`lvelasymax' `lvelasymin')
          legend(pos(${gvpos}) ring(0) col(1) colgap(*1) symxsize(*1) 
               order(1 "me (to % changes in HC rental rates)" 
                     2 "ies (to % changes in HC rental rates)" 
                     3 "me' (to % changes in wages)" 
                     4 "ies' (to % changes in wages)"));
   graph export ${gpicsdir}/eps_elas.eps, replace;

   line me0 me2 t if t<=`lvlastT', graphregion(color(white))
          lpattern("l" "_" "-" "_.")  xtitle("Age") ytitle("Elasticity")
          xscale(range(18 `lvlastT')) xlabel(20(10)60 `lvlastT',grid) 
          yscale(range(`lvelasymin' `lvelasymax')) ylabel(0(`lvelasstep')`lvelasymax' `lvelasymin')
          legend(pos(${gvpos}) ring(0) col(1) colgap(*1) symxsize(*1) 
               order(1 "me (to % changes in HC rental rates)" 
                     2 "me' (to % changes in wages)" ));
   graph export ${gpicsdir}/eps_elas_me.eps, replace;
   line me0 t if t<=`lvlastT', graphregion(color(white))
          lpattern("l" "_" "-" "_.")  xtitle("Age") ytitle("Elasticity") 
          xscale(range(18 `lvlastT')) xlabel(20(10)60 `lvlastT',grid) 
          yscale(range(`lvelasymin' `lvelasymax')) ylabel(0(`lvelasstep')`lvelasymax' `lvelasymin');
   graph export ${gpicsdir}/eps_elas_me0.eps, replace;
   line me2 t if t<=`lvlastT', graphregion(color(white))
          lpattern("l" "_" "-" "_.")  xtitle("Age") ytitle("Elasticity")
          xscale(range(18 `lvlastT')) xlabel(20(10)60 `lvlastT',grid) 
          yscale(range(0 1.5)) 
          ylabel(0(0.3)1.5 0.1 0.2);
          *yscale(range(`lvelasymin' `lvelasymax')) ylabel(0(`lvelasstep')`lvelasymax' `lvelasymin');
   graph export ${gpicsdir}/eps_elas_me2.eps, replace;
   line me2 t if t>=20 & t<=`lvlastT', graphregion(color(white))
          lpattern("l" "_" "-" "_.")  xtitle("Age") ytitle("Elasticity")
          xscale(range(18 `lvlastT')) xlabel(20(10)60 `lvlastT',grid) 
          yscale(range(0 1.5)) 
          ylabel(0(0.3)1.5 0.1 0.2);
          *yscale(range(`lvelasymin' `lvelasymax')) ylabel(0(`lvelasstep')`lvelasymax' `lvelasymin');
   graph export ${gpicsdir}/eps_elas20_me2.eps, replace;
   line me2 t if t>=20 & t<=`lvlastT', graphregion(color(white))
          lpattern("l" "_" "-" "_.")  xtitle("Age") ytitle("Elasticity")
          xscale(range(18 `lvlastT')) xlabel(20(10)60 `lvlastT',grid) 
          yscale(log);
   graph export ${gpicsdir}/eps_elas20_me2_log.eps, replace;


   * output ies numbers;
   keep if mod(t,5)==0;
   sort t;
   *outsheet t me0 ies0 me2 ies2 using ${gpicsdir}/_elas0.csv, c replace;
   export excel t me0 ies0 me2 ies2 using "table4paper${alt_name}", sheet("elas${alt_name}") sheetreplace firstrow(variables);

restore;


/*********************************************************
* calculate total lfpr response to change in pw;
*********************************************************/

local lvmoments_elast lfpr lnw lnw_fd `lvAH';
local lvmoments_elast_b "";
foreach iv of local lvmoments_elast {;
   local lvb2 `lvmoments_elast_b' `iv'_b;
   local lvmoments_elast_b `lvb2';
};

preserve;
   foreach iv of local lvmoments_elast {;
      replace `iv' = `iv' - `iv'_b;
   };
   tab t;
   keep `lvmoments_elast' `lvmoments_elast_b' elast;
   collapse (sum) `lvmoments_elast' `lvmoments_elast_b', by(elast);
   sort elast;
   foreach iv of local lvmoments_elast {;
      gen pd_`iv' = `iv' / `iv'_b;
   };
   sort elast;
   compress;
   save ${gdtadir}/dta_MSM_elas_overall.dta, replace;
restore;

preserve;
   gen byte elast_pos = .;
   replace elast_pos = 0 if t<elast;
   replace elast_pos = 1 if t>elast;
   replace elast_pos = 2 if t==elast;

   foreach iv of local lvmoments_elast {;
      replace `iv' = `iv' - `iv'_b;
   };
   keep `lvmoments_elast' `lvmoments_elast_b' elast elast_pos;
   collapse (sum) `lvmoments_elast' `lvmoments_elast_b', by(elast elast_pos);
   reshape wide `lvmoments_elast' `lvmoments_elast_b', i(elast) j(elast_pos);
   sort elast;
   foreach iv of local lvmoments_elast {;
      foreach jv in 0 1 2 {;
         gen pd_`iv'`jv' = `iv'`jv' / `iv'_b`jv';
      };
   };
   sort elast;
   merge 1:1 elast using ${gdtadir}/dta_MSM_elas_overall.dta, nogenerate;
   sort elast;
   compress;
   save ${gdtadir}/dta_MSM_elas_overall.dta, replace;

   * plot graphs now;
   *foreach iv of local lvmoments_elast {;
   foreach iv in lfpr {;
      if ("`iv'"=="lfpr") {;
         local lvytitle "Labor Force Participation Rate";
         local lvpos1 = 4;
         local lvpos2 = 4;
      };
      else if ("`iv'"=="lnw") {;
         local lvytitle "Log Wages";
         local lvpos1 = 10;
         local lvpos2 = 6;
      };
      else if ("`iv'"=="lnw_fd") {;
         local lvytitle "Log Wages (FE)";
         local lvpos1 = 10;
         local lvpos2 = 6;
      };
      else if ("`iv'"=="C") {;
         local lvytitle "Consumption";
         local lvpos1 = 10;
         local lvpos2 = 6;
      };
      else {;
         local lvytitle "`iv'";
         local lvpos1 = 4;
         local lvpos2 = 4;
      };

      line `iv' `iv'0 `iv'1 `iv'2 elast if elast<=${gvt_elas}, lpattern("l" "_" "-" "-.")  xtitle("Shock age") ytitle("`lvytitle'") graphregion(color(white))
          xlabel(20(10)${gvt_elas} ${gvt_elas},grid)
          legend(pos(`lvpos1') ring(0) col(1) colgap(*1) symxsize(*1) order(1 "Overall" 2 "Before t" 3 "After t" 4 "At t"));
      graph export ${gpicsdir}/eps_elas_response_`iv'.eps, replace;

   };
restore;


cap log close;

